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SECTION  1 


INTRODUCTION 


Seismic  calibration  and  yield  estimation  for  the  Novaya  Zemlya  tost  site 
(NZ)  are  still  matters  of  some  interest  and  continued  research;  if  underground  nuclear 
testing  resumes  in  the  former  Soviet  Union,  it  is  likely  to  be  at  the  northern  Novaya 
Zemlya  site  (NNZ).  (No  tests  have  been  performed  at  the  southern  site  (SNZ)  since 
1975,  e.g..  Lay,  1991.)  Accurate  magnitude-log  yield  calibration  for  NNZ  has  been 
hindered  in  the  past  by  the  lack  of  yield  information  to  accompany  a  wide  variety  of 
published  seismic  magnitude  data.  Recently,  yield  data  for  42  underground  nuclear 
tests  conducted  at  NZ  have  been  provided  to  DARPA  by  an  official  of  the  former 
Soviet  Union.  The  focus  of  this  report  is  on  comparing  these  yields  with  previous  esti¬ 
mates  and  on  computing  estimates  of  the  multivariate  magnitude-log  yield  calibration 
parameters. 

Teleseismic  body  wave  magnitudes,  mj,  have  been  given  by  the  International 
Seismic  Centre  (ISC),  Sykes  and  Ruggi  (1988),  Burger  et  al.  (1986),  and  Chan  et 
al.  (1988b).  Nuttli  (1988)  and  Sykes  and  Ruggi  (1988)  provide  sets  of  m\,[Lg)  and  Ms 
magnitudes,  respectively.  Chan  et  al.  (1988a,b)  and  Lilwall  and  Marshall  (1986)  give 
magnitudes  based  on  P  and  P’P’.  These  data  and  yield  estimates  based  on  these  data 
and  others  are  reviewed  by  Lay  (1991).  More  recently,  Ringdal  and  Fyen  (1991)  have 
published  world-wide  m*  values,  RMS  Lg  and  P  coda  values  recorded  at  NORSAR,  and 
RMS  Lg  values  recorded  at  Grafenberg  for  18  events  at  NNZ  between  29  September 
1976  and  24  October  1990.  Israelsson  (1992)  has  published  a  set  of  Soviet  network 
magnitudes  based  on  RMS  amplitudes  in  time  windows  corresponding  to  initial  P,  P 
coda,  S  with  coda,  and  Lg  for  17  of  the  same  18  events. 

Epicentral  locations  of  many  of  the  tests  were  determined  by  Lilwall  and 
Marshall  (1986).  Others  were  determined  at  the  Center  for  Seismic  Studies  (CSS) 
and  the  ISC.  Israelsson  (1992)  provides  a  review  of  the  epicentral  locations  of  21  of 
the  tests  at  NNZ.  Due  to  the  complicated  topography  of  the  northern  site,  these  tests 
were  mainly  emplaced  in  near-horizontal  tunnels  in  the  mountains.  Leith  et  al.  (1990) 
and  Lay  (1991)  provide  reviews  of  the  test  site  geology,  topography,  tectonic  release 
and  complex  surface  interactions  and  propagation  effects. 

The  Soviet  yield  data  used  in  our  study  was  presented  by  Victor  Mikhailov 
of  the  former  Soviet  Union  at  a  meeting  in  Norway  in  September  1991.  Several 
apparent  discrepancies  in  the  data  were  pointed  out  to  us  by  Richards  (1992).  More 
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recent  discussions  with  Adughkin  et  al.  (1992)  at  the  14th  Annual  PL/DARPA  Seismic 
Research  Symposium  in  Tucson  have  resolved  some  of  the  discrepancies.  However, 
they  brought  into  question  the  accuracy  of  the  data.  Thus,  the  magnitude-log  yield 
results  presented  here  are  preliminary  and  contingent  on  the  accuracy  of  the  yield 
data. 


The  remainder  of  this  report  is  organized  as  follows.  In  Section  2  we  present 
the  Soviet  yield  data  and  compare  it  to  previous  yield  estimates  based  on  observed 
seismic  magnitudes  and  magnitude-log  yield  relations  transported  from  other  sites. 
We  discuss  the  discrepancies  associated  with  the  yield  data.  In  Section  2  we  also 
provide  the  magnitude  data,  used  in  our  analysis  (from  Ringdal  and  Fyen,  1991,  and 
Israelsson,  1992). 

In  Section  3  we  provide  the  statistical  background  needed  for  this  study.  In 
Section  3.1  we  present  multivariate  regression  analysis  and  in  Section  3.2  we  present  a 
95%  confidence  interval  for  the  yield  of  a  new  event.  The  confidence  interval  is  based 
on  multivariate  data  and  the  treatment  of  all  model  parameters  as  unknown.  It  was 
originally  derived  by  Brown  (1982)  and  has  been  used  previously  by  Shumway  and  Der 
(1990)  to  compute  yield  intervals  for  Semipalatinsk  explosions.  We  also  describe  how 
TTBT  compliance  may  be  tested  and  define  the  F-number  in  terms  of  the  confidence 
interval.  In  Section  3.3  we  present  an  outlier  test,  which  we  later  apply  to  the  data. 

In  Section  4  we  provide  the  results  of  calibration,  outlier,  and  yield  estima¬ 
tion  for  the  magnitude  data  from  Ringdal  and  Fyen  (1991).  In  Section  4.1  we  provide 
estimates  of  the  calibration  parameters,  including  the  random  error  correlation  coef¬ 
ficients.  In  Section  4.2  we  discuss  which  events  have  residuals  that  are  identified  as 
outliers.  In  Section  4.3  we  provide  estimates  of  the  F-number  means  and  standard 
deviations  for  various  magnitude  combinations,  obtained  by  jackknifing.  In  Section  5 
we  repeat  this  analysis  for  the  Soviet  network  magnitudes  for  Israelsson  (1992).  Last, 
in  Section  6  we  provide  some  concluding  remarks  and  a  review  the  best  calibration 
parameter  estimates. 
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SECTION  2 


DATA 


2.1  Novaya  Zemlya  Yield  Data. 


Within  the  last  year,  yield  data  for  42  underground  nuclear  tests  conducted 
at  NZ  were  provided  to  DARPA  by  an  official  of  the  former  Soviet  Union.  The 
original  bar  graph  of  the  yield  data  is  shown  in  Figure  1.  The  year  and  number  of 
events  tested  in  that  year  (in  parentheses)  are  labeled  along  the  i-axis,  and  the  month 
and  accumulated  number  of  events  tested  during  a  particular  month  are  labeled  along 
the  left  and  right  edges  of  the  j/-axis,  respectively.  The  dates  of  each  event  are  listed 
at  the  top  of  the  graph.  We  have  determined  the  yields  by  digitizing  the  bar  lengths 
and  comparing  them  to  the  length  of  the  150  KT  scale  at  the  left  edge  of  the  graph. 

Table  1  lists  the  dates  provided  with  the  original  data,  our  determination 
of  the  yields,  Y(FGM),  as  well  as  previous  estimates  of  some  of  the  yields  by  Nuttli 
(1988),  Y(Nuttli),  Sykes  and  Ruggi  (1988),  Y(SR),  and  Burger  et  al.  (1986),  Y(BBL). 
The  yield  estimates  by  Nuttli  (1988)  were  based  on  a  quadratic  fit  of  Lg  magnitudes  to 
the  log  yields  of  NTS  events.  The  fit  was  then  transported  to  NZ  with  no  corrections. 
The  yields  from  Sykes  and  Ruggi  (1988)  were  estimated  from  a  linear  regression  model 
for  mj,  assuming  bias  corrections  relative  to  NTS  and  Amchitka.  The  yields  from 
Burger  et  al.  (1986)  were  estimated  by  scaling  relative  to  an  Amchitka  event.  These 
and  other  yield  estimates  of  NZ  events  have  been  compiled  previously  by  Lay  (1991); 
we  present  them  here  for  comparison.  Jih  and  Wagner  (1992)  have  also  estimated  the 
yields  of  28  of  these  events  based  on  path  corrected  short-period  teleseismic  P-wave 
amplitudes  and  experience  at  Semipalatinsk  and  NTS. 

There  are  several  points  regarding  the  new  yield  set  that  should  be  noted. 
First,  there  are  breaks  in  7  of  the  8  longest  bars.  These  events  are  larger  that  1  MT, 
as  confirmed  by  the  yield  estimates,  but  the  bars  do  not  provide  accurate  measures  of 
the  yields.  Second,  our  estimates  of  the  very  small  events  (12,  20,  21,  37)  have  limited 
accuracy;  in  fact,  we  do  not  currently  know  the  accuracy  with  which  any  of  the  yields 
were  recorded  on  this  graph.  Third,  the  graph  in  Figure  1  shows  a  small  event  on  27 
July  1972,  which  is  not  present  in  any  of  the  other  data  sets  we  have  examined.  We 
have  learned  from  Richards  (1992)  that  this  test  may  have  never  occurred.  Adushkin 
et  al.  (1992)  have  confirmed  that  this  test  did  not  occur. 
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Figure  1.  Original  bar  graph  of  Soviet  yield  data  for  42  underground  nuclear 
tests  conducted  at  Novaya  Zemlya. 

Fourth,  the  two  events  on  18  October  1975  are  known  to  be  double  explosions 
at  SNZ  (Hurley,  1977;  Subhash  and  Choudhury,  1979;  Burger  et  al.,  1986,  Chan  et  al., 
1988a),  and  all  available  yield  estimates  suggest  that  these  were  large  events,  on  the 
order  of  a  megaton  or  more.  These  events  are,  however,  characterized  as  very  small 
on  the  bar  graph  (Figure  1),  This  discrepancy  was  first  pointed  out  to  us  by  Richards 
(1992).  Adushkin  et  al.  (1992)  have  confirmed  that  the  yields  were,  in  fact,  on  the 
order  of  a  megaton  or  more. 

Fifth,  Lilwall  and  Marshall  (1986)  and  Stewart  and  Marshall  (1988)  sug¬ 
gested  that  the  event  on  11  October  1980  was  also  a  double  explosion,  with  shots 
separated  by  approximately  7  krn.  They  determined  that  the  ratio  of  the  yields  was 
roughly  a  factor  of  0.35.  This  is  consistent  with  our  yield  determinations  given  in 
Table  1,  Y(FGM),  for  these  events. 
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Table  1,  Yield  determinations  and  estimates  of  Novaya  Zemlya  tests. 


Event 

Date 

Y  (FGM)  i 

Y  (Nuttli)  i 

I  Y  (SR) 

Y  (BBL) 

1  j 

640918 

29 

1  2.51 

1  21 

2I 

641025 

40 

16. 4' 

8 

1 

3\ 

661027 

603 

!  644 

422 

!  600 

4i 

661027 

485 

i  _ 

1 

5 

671021 

184 

180 

93 

61 

6 

126 

- 1 

7 

681107 

334 

253 

1  19 

110 

8 

691014 

302 

399 

140 

183 

9 

691014 

219 

1  0 

701014 

>1000 

1970 

1001 

1714 

1  1 

710927 

>1000 

586 

973 

1  2 

720727 

14 

13 

720828 

1183 

580 

329 

426 

1  4 

730912 

>1000 

3510 

2099 

2824 

1  5 

730927 

102 

129 

100 

36 

1  6 

731027 

>1000 

4990 

4055 

3886 

1  7 

740829 

>1000 

1110 

497 

629 

1  8 

741102 

>1000 

2840 

1624 

19 

^  750823 

152 

690 

477 

604 

20 

751018 

10 

2220 

1281 

1166 

21 

751018 

10 

22 

751021 

>1000 

600 

497 

554 

23 

760929 

131 

91 

70 

1  _ 

24 

761020 

22 

19 

13 

25 

770901 

134 

122 

55 

26 

771009 

27 

10 

4 

780810 

121 

91 

89 

28 

780927 

125 

61 

44 

29 

790924 

144 

81 

55 

30 

791018 

120 

79 

70 

31 

801011 

139 

76 

55 

32 

801011 

30 

33 

811001 

136 

116 

113 

34 

821011 

107 

79 

44 

35| 

830818 

105 

145 

89 

830925 

97 

99 

70 

Bn 

840826 

10 

38 

841025 

101 

_ 1 

89 

39 

870802 

121 

70 

40 

880507 

97 

41 

881204 

112 

- 1 

42 

901024 

64 

5 


Last,  many  of  the  yields  are  consistent  with  previous  estimates;  however, 
there  are  several  events  for  which  the  yields  differ  significantly.  Adushkin  et  al.  (1992) 
suggested  that  the  accuracy  of  the  bar  graph  in  Figure  1  is  questionable.  Clearly,  the 
accuracy  of  these  yield  data  needs  to  be  resolved,  and  we  are  pursuing  this  matter 
currently.  Hence,  although  there  are  no  known  discrepancies  associated  with  the 
events  used  to  obtain  the  calibration  results,  given  below,  our  results  are  preliminary 
and  contingent  on  the  accuracy  of  the  original  yield  data. 


2.2  Magnitude  Data. 

Ringdal  and  Fyen  (1991)  have  published  a  set  of  magnitudes,  based  on 
NORSAR  Lg  and  P  coda,  Grafenberg  Lg,  and  a  world-wide  mj,  for  18  undergroimd 
nuclear  tests  conducted  between  29  September  1976  and  24  October  1990  at  NNZ. 
(See  Ringdal  and  Fyen,  1991,  and  references  therein  for  descriptions  of  the  arrays  and 
the  methods  used  to  compute  the  magnitudes.)  Table  2  lists  the  event  numbers,  dates, 
yields  and  magnitudes  of  those  tests.  Low  SNR  at  Grafenberg  for  event  26  did  not 
allow  for  a  reliable  RMS  Lg  determination.  Also,  no  NORSAR  data  are  available  for 
event  28,  and  event  37  was  too  small  to  provide  reliable  magnitude  determinations  from 
either  array.  The  magnitudes  corresponding  to  the  double  explosion  on  11  October 
1980  (events  31  and  32)  may  be  biased  from  interference  effects.  We  have  associated 
the  magnitudes  with  the  larger  of  the  two  yields;  using  the  root  sum  of  the  squares  of 
the  combined  yields  (142  KT)  produces  an  insignificant  change  in  the  results. 

Israelsson  (1992)  has  published  a  set  of  Soviet  network  magnitudes,  which 
we  will  also  consider  in  our  study.  The  network  is  comprised  of  nine  stations,  Ap¬ 
atity  (APA),  Arti  (ARU),  Bodaybo  (BOD),  Chusal  (CHS),  Norilsk  (NRI),  Novosibirsk 
(NVS),  Obninsk  (OBN),  Talaya  (TLY),  and  Uzhgorod  (UZH).  (See  Israelsson,  1992, 
for  the  locations  of  these  stations,  a  description  of  the  seismometers,  and  further  de¬ 
tails  on  the  results.)  Analog  recordings  of  111  waveforms  for  21  explosions  from  the 
nine  stations  were  hand-digitized.  Although  his  study  emphasized  RMS  Lg  in  the 
group  velocity  window  between  3. 1-3.7  km/s,  he  also  computed  RMS  amplitudes  cor¬ 
responding  to  initial  P  in  the  time  window  from  onset  to  20  s  after  onset,  P  coda  in 
the  time  window  20  s  after  onset  to  15  s  before  the  expected  S  arrival,  and  S  with 
coda  in  the  time  window  15  s  prior  to  the  expected  S  arrival  to  a  group  velocity  of 
3.7  km/s. 


RMS  amplitudes,  computed  from  bandpass  filtered  traces,  were  corrected 
for  noise  and  normalized  to  the  OBN  instrument.  Station  magnitudes  were  obtained 
by  taking  the  logarithm,  and  then  shifted  by  a  constant,  for  each  station  and  phase. 
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Table  2.  Magnitude  data  from  Ringdal  and  Fyen  (1991)  for  18  tests  at  NNZ 


Event 

Date 

Y  (kt) 

NAO  Lg 

GRFLg  ■ 

NAO  P  coda  [ 

world-wide  mb 

23 

760929 

131 

5.770 

5.7991 

5.732' 

5.77 

24 

761020 

22 

5.071 

5.0221 

4.969 

4.89 

25 

770901 

134 

5.757 

5.8721 

5.750 

5.71 

26 

771009 

27 

4.845 

0.0001 

4.637; 

4.51 

27 

780810 

im 

5.783 

5.759i 

5.9521 

6.04 

28 

780927 

125 

0.000 

5.6601 

O.OOOi 

5.68 

29 

790924 

144 

5.779 

5.8251 

5.782 

5.80 

30 

791018 

120 

5.737 

5.664| 

5.821  i 

5.85 

31 

801011 

139 

5.784 

5.7321 

5.7761 

5.80 

33 

811001 

136 

5.782 

5.783i 

5.8821 

5.91 

34 

821011 

107 

5.603 

5.5851 

5.551! 

5.52 

35 

830818 

105 

5.807 

5.7391 

5.7691 

5.84 

36 

830925 

1  97 

5.797 

5.7771 

5.723! 

5.71 

38 

841025 

1011  5.805 

5.8371 

5.7431 

5.77 

39 

870802 

121 

5.806 

S.810I 

5.7691 

5.71 

40 

880507 

97 

5.7191  5.654i 

5.6141 

5.52 

41 

881204 

112 

5.800 

5.811; 

5.822i 

42 

901024 

64 

5.6051  5.550; 

5.618! 

5.60 

Table  3.  Magnitude  data  from  Israelsson  (1992)  for  17  tests  at  NNZ. 


Event 

Date 

Y  (kt) 

mb(P  coda) 

mb(S  coda) 

23 

7609291 

131 

5.780 

5.695; 

5.651 

5.705 

24 

761020 

22 

5.048 

5.031 

4.993 

25 

770901 

134 

5.803 

5.821! 

5.750 

5.745 

26 

771009 

27 

4.788 

4.819' 

4.698 

4.735 

27 

780810 

121 

5.820 

5.825! 

5.796 

5.786 

28 

780927 

125 

5.640 

5.568; 

5.616 

5.597 

29 

790924 

144 

5.863 

5. 8411 

5.745 

5.746 

30 

791018 

120 

5.813 

5.7561 

5.725 

5.720 

31 

801011 

139 

5.764 

5.749i 

5.803 

5.750 

33 

811001 

136 

5.853 

5.8271 

5.798 

5.798 

34 

821011 

107 

5.586 

5.587i 

5.581 

35 

830818 

105 

5.792 

5.779i 

5.813 

36 

830925 

97 

5.728 

5.8111 

5.793 

5.809j 

38 

841025 

101 

5.713 

5.700i 

5.678 

39 

870802 

121 

5.835 

5.89li 

5.827!  5.797 

40 

880507 

97 

5.644 

5.716!  5.714 

41 

881204 

5.7741  5.738! 

5.7691  5.810 

to  normalize  them  relative  to  the  mean  NAO  Lg  magnitude,  based  on  the  events 
with  magnitude  between  5.603-5.807  in  Table  2.  Station  magnitudes  are  assumed 
to  be  a  linear  combination  of  the  network  magnitude,  station  correction,  and  zero 
mean  Gaussian  error  terms.  Network  magnitudes  were  obtained  by  solving  the  over- 
determined  system  of  equations  for  the  network  magnitude,  station  correction,  and 
standard  deviation  using  a  least  squares  method.  The  resulting  network  magnitudes 
based  on  the  four  RMS  amplitudes  are  given  in  Table  3. 

This  set  of  17  events  includes  all  but  event  42  on  24  October  1990  of  the  set 
from  Ringdal  and  Fyen  (1991).  The  17  Soviet  network  magnitudes  are  based  on  mea¬ 
surements  varying  from  three  to  eight  stations.  There  were  four  other  events  analyzed 
by  Israelsson,  events  2,  4,  7,  and  37,  but  no  network  magnitudes  were  computed  for 
the  first  three  because  of  uncertainties  in  the  instrument  characteristics  at  the  time 
of  these  events.  Network  magnitudes  were  computed  for  event  37  on  26  August  1984; 
however,  they  are  based  on  measurements  from,  at  most,  two  stations  and  at  least  one 
of  the  measurements  was  characterized  as  having  poor  SNR.  Note  also  that  the  yield 
for  this  event  is  represented  by  a  dot  in  Figure  1,  which  leads  to  considerable  uncer¬ 
tainty  in  the  determined  yield.  Thus,  we  have  omitted  this  event  from  the  analysis 
below  because  neither  the  yield  nor  the  magnitudes  for  this  event  are  very  reliable. 


8 


SECTION  3 


STATISTICAL  BACKGROUND 


3.1  Calibration. 


In  this  section,  we  present  the  standard  linear  magnitude-log  yield  model, 
for  a  vector  of  magnitude  observations,  and  describe  how  the  model  parameters  may 
be  estimated  from  data.  Let  m  =  (mi,m2, . . .  ,mp)'  denote  the  vector  of  p  observed 
seismic  magnitudes,  where  the  prime  superscript  denotes  the  vector  transpose,  and 
w  denote  the  log  yield  of  an  event.  The  linear  magnitude-log  yield  regression  model 
may  then  be  expressed  as 


m  =  a  f  bu;  +  e 


E{e)  =  0  (1) 

Cow(e)  =  E, 

where  the  p  x  1  vectors,  a,  b  and  e,  denote  the  intercept  and  slope  parameters,  and 
the  vector  of  random  seismic  errors,  respectively.  We  will  assume  that  e  has  a  normal 
distribution  with  zero  mean  and  covariance  matrix  denoted  by  E.  It  is  implicitly 
assumed  in  this  model  that  the  log  yields  are  fixed  observations  without  error. 

Direct  estimates  of  the  calibration  parameters,  a,  b  and  E,  may  be  com¬ 
puted  if  a  set  of  measured  seismic  magnitudes  and  corresponding  log  yields,  m,  and 
Wi  (for  i  =  l,...,n.  events),  are  available  from  historical  underground  nuclear  tests. 
For  example,  the  vector  mj  can  represent  the  magnitudes  based  on  NAO  Lg,  GRF  Lg, 
NAO  P  coda,  and  world-wide  mj  for  the  i  =  1, . . . ,  18  events  listed  in  Table  2.  Given 
n  previously  observed  values  of  m,,  for  which  the  Wi  are  known,  unbiased  estimates 
of  b,  a  and  E  are  given  by 


n 

-  rn)(u;<  -  «J) 


•=i 


(2) 
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A 

a  =  m  —  hw 
1  n 

±  =  — -a-bu;i)(in< -a-bty,)', 

”  ^  t=i 

where  in  and  w  are  the  usual  sample  means,  given  by 


(3) 

(4) 


1  " 

".=1 

(5) 

1  " 

(6) 

(See,  e.g.,  Fuller,  1987.)  The  estimates  of  the  variances  and  correlation  coefficients,  in 
terms  of  the  estimates  of  the  covariance  matrix  elements,  are  given  by 


=  (S]«;  i  = 


Pjk 


E, 


it _ . 


v/e,,.e 


;  j,k  =  l,...,p. 


(7) 

(8) 


kk 


3.2  Yield  Estimation. 


Here  we  follow  the  analysis  of  Shumway  and  Der  (1990)  to  compute  a  clas¬ 
sical  confidence  interval  for  the  yield  of  a  new  event.  To  motivate  the  calculation  of 
the  confidence  interval,  we  first  consider  the  minimum  variance,  unbiased  estimator 
of  the  log  yield 


b'E  ‘(m-a) 

tw  = - -  — - 

O  d 

where  m  is  the  magnitude  vector  measured  for  the  new  event.  Assuming  that  m 
has  a  multivariate  normal  distribution,  a  confidence  interval  for  the  log  yield  may  be 
computed  and  a  test  of  hypothesis  that  the  yield  is  below  an  arbitrary  threshold  may 
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be  established.  Unfortunately,  the  calibration  parameters  are  unknown,  rendering 
this  analysis  infeasible.  It  is  tempting  to  replace  the  unknown  parameters  with  their 
estimators,  a,  b  and  E  in  equation  (9);  however,  the  resulting  distribution  of  the  log 
yield  estimator  is  not  well  known,  and  the  confidence  interval  cannot  be  calculated 
readily. 

An  alternative  confidence  interval,  considered  by  Shumway  and  Der  (1990), 
was  originally  derived  by  Brown  (1982),  and  has  been  discussed  by  Anderson  (1984), 
Oman  (1988)  and  others.  Before  deriving  the  confidence  interval  it  is  useful  to  express 
the  multivariate  model  for  n  observations  as 


M  =  B'W  +  E,  (10) 

where  M  =  (mj, . . . ,  m„)  and  E  =  (ci, . . .  ,e„)  are  p  x  n  vectors,  B'  —  (a,  b)  is  a  p  x  2 
vector,  and 


The  minimum  variance,  unbiased  estimator  of  B  is 

B  =  C"^WM', 


(11) 


(12) 


where 


c  =  ww'  = 


n  nw  \ 
nuJ  y 


(13) 


It  is  straightforward  to  show  that  the  expressions  in  (2)  and  (3)  may  be  recovered 
from  (12).  Also,  the  unbiased  estimator  of  E  is  given  by 


E  =  (n-2)-‘ (M-BW)(M-BW)'.  (14) 

Using  the  assumptions  that  e  ~  iV(0,E)  and  that  the  random  error  vectors  are  inde¬ 
pendent  for  different  events,  it  may  be  shown  that  B  ~  A'^(B,E  (gi  C“‘),  i.e.,  has  a 
normal  distribution  with  mean  B  and  covariance  matrix  E0C~*,  where  ®  denotes  the 
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Koenecker  product.  Using  Theorems  4.3.3  and  7.2.2  of  Anderson  (1984,  pp.  130,  249), 
it  may  also  be  shown  that  (n  —  2)X)  ~  W(E,n  —  2),  i.e.,  has  a  Wishart  distribution 
with  n  —  2  degrees  of  freedom. 

To  compute  a  confidence  interval,  note  that  for  a  new  observed  magnitude 
vector  m,  associated  with  unknown  log  yield  w,  the  multivariate  residual  is  distributed 
as 


m  —  a  —  biy  ~  (1  +  9(u;))E),  (15) 

where  q{w)  =  (1, u;)C”‘(l, u;)'.  Now  let 


=  (m  -  a  -  bio)'  [(l  +  9(io))Ej  (m-a-bto).  (16) 

Using  Theorem  5.2.2  of  Anderson  (1984,  p.  163),  it  may  be  shown  that 


2  n-p- 1 
p(n-2) 


(17) 


where  Fp^n-p-i  ‘s  the  F  distribution  with  p  and  n  —  p  —  I  degrees  of  freedom.  (The 
distribution  of  T*  is  referred  to  as  Hotellings  T*  distribution;  T*  is  the  multivariate 
analog  of  the  square  of  a  random  variate  t  which  has  Student’s  t  distribution.) 

A  100(1  -  a)%  confidence  interval  may  be  obtained  for  the  log  yield  lo  by 
inverting  the  quadratic  inequality 


r'  <  (18) 

where  Fp,„_p_i(a)  is  the  (1  -  a)  percentile  of  the  F  distribution.  Denoting  the  matrix 
elements  of  C'*  by  and 

C  =  h'±-%  -  C^Xn-2H  (19) 

d  =  b'E-‘(ni-a)  +  c‘*T,V2(«)  (20) 

e  =  (m-a)'E"‘(m-a)  -  (1 +  c")Tp*„_2(a),  (21) 
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the  endpoints  of  the  interval  are  given  by 


d  ^  yjd}  —  ct 
c  c 

Brown  (1982)  has  noted  that  the  interval  may  be  open  and  is  not  necessarily  real. 
Open  intervals  occur  when  the  quadratic,  T*,  is  concave  down  as  a  function  of  w.  The 
roots  of  the  quadratic  inequality  for  w  may  also  be  complex.  These  situations  occur  if 
the  new  magnitude  information  sufficiently  contradicts  the  calibration  data,  or  if  the 
spread  in  the  calibration  data  is  too  small.  Examples  of  these  pathologies  will  be  seen 
in  the  results  below. 

Once  the  confidence  interval  has  been  computed,  the  midpoint  of  the  interval 
is  commonly  used  as  a  log  yield  point  estimate.  A  yield  confidence  interval  and  point 
estimate  may  be  obtained  by  exponentiating  the  log  yield  interval  endpoints  and 
midpoint  to  the  tenth  power.  The  yield  and  log  yield  point  estimates  are  likely  to  be 
biased,  i.e.,  E\dlc\  ^  w  and  ^  y  =  lO",  in  general. 

A  test  of  the  hypothesis  that  the  yield  is  in  compliance  with  a  treaty  thresh¬ 
old  may  be  established  of  the  form:  Reject  the  hypothesis  if  the  lower  endpoint  of  the 
confidence  interval  for  the  yield  is  greater  than  the  treaty  threshold.  If  the  central 
value  of  the  100(1  -  a)%  confidence  interval  is  unbiased,  the  significance  level.  A,  of 
this  test  would  be  a/2.  In  general,  however,  the  central  value  is  a  biased  estimator  of 
the  yield.  Thus,  it  can  only  be  stated  that  0  <  A  <  a.  (If  we  had  an  unbiased  estimate 
of  the  log  yield,  a  strict  a/2  significance  level  test  could  be  established,  although  it 
could  not  be  expressed  in  closed  form  unless  the  distribution  was  also  known.  Treating 
the  slope  parameters  as  unknown  is  the  fundamental  cause  of  both  complications.) 

We  define  the  F-number  here  such  that  the  yield  interval  may  be  expressed 
as  [Y / F,  y  X  F),  where  V  is  the  central  value  of  the  interval.  Thus,  it  is  given  by 

F  =  (23) 

Nicholson  et  al.  (1991)  also  define  an  F-number  in  terms  of  a  confidence  interval.  The 
confidence  interval  they  derive  is  based  on  a  Bayesian  approach,  which  treats  the  slope 
and  intercept  parameters  as  unknown,  while  treating  the  covariance  matrix  as  known. 

Note  that  this  definition  is  intrinsically  different  from  the  one  presented 
by  Alewine  et  al.  (1988).  The  F-number  given  here  is  a  random  variable,  while  the 
F-number  defined  by  Alewine  et  al.  (1988)  is  a  fixed  quantity,  defined  to  be  the 
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value  of  the  actual  yield,  such  that  the  probability  of  detecting  a  treaty  violation 
is  0.5,  divided  by  the  treaty  yield  threshold.  The  definition  of  the  F-number  as  a 
fixed  number,  given  by  Alewine  et  al.  (1988),  has  an  appealing  interpretation  and, 
in  principle,  could  be  applied  here.  Unfortunately,  determination  of  this  F-number 
depends  on  the  distribution  of  the  log  yield  estimator  and,  implicitly,  on  the  unknown 
slope,  intercept  and  covariance  parameters.  Thus,  we  will  adopt  the  definition  given 
here  and,  later,  estimate  the  mean  and  standard  deviation  of  the  random  variable,  F, 
by  applying  a  jackknife  procedure  to  the  data. 


3.3  Outlier  Detection. 

The  statistic  presented  in  the  previous  section  can  also  be  used  to  de¬ 
termine  if  a  particular  event  in  the  calibration  data  is  an  outlier,  i.e.,  to  test  whether 
the  residual  of  an  event  is  anomalously  large.  To  test  whether  a  given  event  is  an 
outlier,  we  first  perform  the  calibration  with  the  remaining  data  to  obtain  estimates 
of  the  intercepts,  slopes  and  covariance  matrix.  These  estimates  and  the  magnitude 
and  log  yield  of  the  event  in  question  are  then  inserted  into  the  expression  for  T*  in 
(16),  whose  value  we  denote  by  Tg .  Since  we  know  the  distribution  of  T*,  we  can 
compute  the  probability  of  obtaining  a  value  of  T*  greater  than  or  equal  to  Tq  and, 
hence,  establish  a  test  of  hypothesis. 

For  p  =  1  (i.e.,  for  a  one  dimensional  magnitude  vector),  T*  is  proportional 
to  the  ratio  of  the  residual  squared  to  the  sample  random  error  variance,  emd  the 
distribution  reduces  to  The  F  statistic  is  commonly  used  to  test  whether  two 

samples  have  the  same  population  mean.  That  is,  in  fact,  what  we  are  testing,  i.e.,  that 
the  population  mean,  pg,  of  the  residual  m  —  a  —  bu;  for  the  event  in  question  is  equal 
to  the  population  mean,  p,  of  the  remaining  residuals  m|  —  a  —  bu;,,  i  =  1, . . . ,  n  —  1. 
As  the  residual  m  —  a  -  bu;  becomes  large,  the  probability  of  obtaining  a  value  of 
T*  greater  than  or  equal  to  Tq  becomes  small.  For  some  critical  value,  we  reject  the 
hypothesis  that  the  residual  of  the  event  in  question  has  the  same  population  mean 
as  the  remaining  sample  and  call  the  event  an  outlier. 

More  formally,  let  Hq  :  hq  =  fi  he  the  hypothesis  that  the  population  means 
are  equal.  The  null  hypothesis  Hq  is  rejected,  with  significance  level  a,  if  P\T^  > 
Tq]  <  a  or,  equivalently,  if  Tq  >  Tpn_2{a).  This  test  is  equivalent  to  rejecting  Hq  if 
the  log  yield  of  the  event  in  question  is  not  included  in  the  100(1  -  a)%  confidence 
interval  defined  above.  (There  are  other  rigorous  tests  for  outliers  that  could  also  be 
used,  e.g.,  a  likelihood  ratio  criterion.)  We  will  apply  the  test  given  here  to  the  data 
in  the  following  sections. 
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SECTION  4 


RESULTS:  RINGDAL  AND  FYEN  MAGNITUDE  DATA 


4.1  Calibration  Results. 


Using  the  data  in  Table  2,  we  have  computed  dy,  bj,  Oj,  and  pjk  {j,k  - 
1, . . .  ,4).  Since  some  magnitudes  were  not  observed  by  all  of  the  stations,  estimates 
of  the  intercepts,  slopes  and  standard  deviations  were  computed  using  all  available 
magnitude  data  recorded  by  the  particular  station,  while  random  error  correlation 
coefficient  estimates  were  computed  using  only  the  intersection  of  events  that  were 
recorded  in  common  by  both  relevant  stations  or  networks.  That  is,  the  intercept, 
slope  and  standard  deviation  estimates  used  in  the  correlation  estimate  calculation 
were  recomputed  using  only  the  events  for  which  there  were  magnitude  measurements 
from  both  stations.  This  correlation  coefficient  estimator  is  more  robust  than  one  that 
uses  intercept,  slope  and  standard  deviation  estimates  computed  from  all  available 
data. 


Figure  2  shows  scatter  plots  of  the  data  and  lines  of  best  fit.  Three  lines  of 
best  fit  were  computed  using  all  available  data  (solid),  all  except  event  26  (dashed), 
and  all  except  events  24  and  26  (dotted).  Events  24  and  26  are  represented  by  the 
solid  square  and  triangle  markers,  respectively.  The  intercept,  slope  and  standard  de¬ 
viation  estimates  are  given  in  the  legends  above  each  frame.  Random  error  correlation 
coefficient  estimates,  using  only  the  events  in  common  between  relevant  magnitude 
types  are  given  in  Tables  4-6  for  these  three  cases.  We  have  presented  these  cases  to 
show  the  effect  of  the  two  small  events  on  the  slope  estimates. 

We  have  several  comments  regarding  these  results.  First,  NAO  and  GRF 
Lg  magnitudes  exhibit  the  least  random  scatter  of  the  four  sets  of  magnitudes.  The 
sample  random  error  standard  deviation  of  NAO  Lg  is  0.101  for  17  events,  but  decreases 
to  0.071  when  event  26  is  omitted,  and  to  0.057  when  both  events  24  and  26  are 
omitted.  For  GRF  Lg,  the  value  is  0.081  for  a  slightly  different  set  of  17  events,  which 
does  not  include  event  26.  It  decreases  to  0.078  when  event  24  is  omitted.  Since  the 
standard  deviations  are  based  on  slightly  different  sets  of  events,  it  is  difficult  to  make 
a  straightforward  comparison.  Thus,  we  have  recomputed  the  results  using  only  the 
16  events  that  are  common  to  all  four  magnitude  types.  Figure  3  shows  scatter  plots 
of  the  data  and  lines  of  best  fit.  The  intercept,  slope  and  standard  deviation  estimates 
are  given  in  the  legends.  For  this  case,  the  sample  random  error  standard  deviation 
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NAO  Pcoda 


Figure  2.  Plots  of  NAO  Lg,  GRF  Lg,  NAO  P  coda,  and  world-wide  mj 
magnitudes  versus  the  yields  of  18  nuclear  explosions  at  NNZ.  The 
solid,  dashed  and  dotted  lines  in  each  frame  represent  the  lines  of 
best  fit  using  all  of  the  available  18  events,  all  except  event  26 
(solid  triangle),  and  all  except  events  26  and  24  (solid  square).  The 
parameter  estimates  are  given  in  the  legends  above  the  plots. 
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Table  4.  Random  error  correlation  coefficient  estimates  based  on  the  magnitude 
data  from  Ringdal  and  Fyen  (1991). 


Magnitude 

NAO  Lg 

GRF  Lg 

NAO  Pcoda 

mj 

NAO  Lg 

1.000 

0.787 

0.844 

0.771 

GRF  Lg 

1.000 

0.441 

0.383 

NAO  Pcoda 

1.000 

0.968 

1.000 

Table  5.  Correlation  coefficient  estimates  excluding  event  26. 


Magnitude 

NAO  Lg 

GRF  Lg 

NAO  Pcoda 

mi, 

NAO  Lg 

1.000 

0.787 

0.664 

0.557 

GRF  Lg 

1.000 

0.441 

0.383 

NAO  Pcoda 

1.000 

0.946 

rrii 

1.000 

Table  6.  Correlation  coefficient  estimates  excluding  event  24  and  26. 


Magnitude 

NAO  Lg 

GRF  Lg 

NAO  Pcoda 

mi. 

NAO  Lg 

1.000 

0.769 

0.579 

0.506 

GRF  Lg 

1.000 

0.346 

0.308 

NAO  Pcoda 

1.000 

0.947 

mi, 

1.000 
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=  3.773 
=  0.959 
=  0.078 


w  (log  kt) 


w  (log  kt) 


Figure  3.  Plots  of  NAO  Lg,  GRF  Lg,  NAO  P  coda,  and  world-wide  m* 
magnitudes  versus  the  yields  of  16  nuclear  explosions  at  NNZ.  The 
same  set  of  events  are  used  for  all  four  magnitude  types.  The 
parameter  estimates  are  given  in  the  legends  above  the  plots. 


for  NAO  Lg  is  slightly  smaller  than  for  GRF  Lg;  both  are  smaller  than  those  for  NAO 
P  coda  and  world-wide  mi,,  the  latter  being  the  largest.  The  correlation  coefficient 
estimates  are  the  same  as  those  obtained  in  Table  5,  with  the  exception  of  the  GRF 
Lg  -mj  correlation  estimate,  which  is  0.325  for  this  case. 

Second,  the  slope  estimates  for  the  four  magnitudes,  based  on  all  of  the 
events  (solid  lines),  are  larger  than  expected,  particularly  for  NAO  P  coda  and  world¬ 
wide  mj  (Figure  2).  As  a  result,  the  intercept  estimates  are  smaller  than  expected. 
Note  that  the  slope  and  intercept  estimates  are  highly  dependent  on  the  two  smallest 
events.  For  example,  when  either  event  26  or  events  24  and  26  are  omitted,  the  slope 
estimates  decrease  significantly  (Figure  2).  This  has  been  referred  to  as  the  “lolly-pop” 
effect  in  regression  analysis  literature.  Sample  random  error  standard  deviations  for  all 
four  magnitudes  also  dccrecise  when  these  events  are  omitted,  and  estimates  of  the  N  AO 
Lg-NAO  P  coda  and  NAO  Lg-mj  correlation  coefficients  change  significantly  (Tables 
4-6).  Since  GRF  Lg  was  not  recorded  for  event  26,  the  corresponding  correlation 
estimates  only  change  when  event  24  is  also  omitted  and,  even  then,  unsubstantially. 

Third,  the  random  error  correlation  coefficient  estimates  are,  in  some  cases, 
remarkably  high.  For  example,  the  NAO  P  coda-m^  correlation  estimates  range  from 
0.946  to  0.968,  depending  on  whether  events  24  and  26  are  included  in  the  analysis 
(Tables  4-6).  Based  on  all  17  events  for  which  there  were  data,  the  NAO  Lg-NAO  P 
coda  and  NAO  Lg-mj  correlation  estimates  are  0.844  and  0.771,  respectively.  These 
values  decrease  when  events  26  and  24  are  omitted,  but  are  still  greater  than  0.5. 
Similarly,  the  GRF  Lg-NAO  Lg  correlation  estimate  is  0.787,  based  on  all  events 
recorded  in  common,  and  0.769  when  event  24  is  excluded.  The  smallest  correlation 
estimates  are  obtained  for  GRF  Lg-NAO  F  coda  and  GRF  Lg-mj. 


4.2  Outlier  Analysis  Results. 


Significant  changes  in  the  calibration  estimates  when  event  26,  in  particular, 
is  excluded  suggests  that  this  event  may  be  an  outlier  in  the  sense  that  the  residu¬ 
als  to  the  lines  of  best  fit  are  anomalously  large.  We  have  applied  the  outlier  test, 
described  in  Section  3.3,  successively  to  all  of  the  events,  and  for  all  possible  mag¬ 
nitude  combinations,  i.e.,  for  all  distinct  combinations  of  p  =  1,2, 3, 4  dimensional 
magnitude  vectors.  For  this  analysis  only  the  events,  for  which  all  magnitudes  of  the 
p  dimensional  vector  were  observed,  were  tested  or  included  in  the  calibration. 

Denoting  the  NAO  Lg,  GRF  Lg,  NAO  P  coda,  and  world-wide  mi,  magnitudes 
by  mi,  m2,  m3,  and  m4,  Table  7  lists  the  probabilities  of  obtaining  values  of  the 


19 


Table  7.  Probability  of  T*  >  Tq  for  events  detected  as  outliers  for  all  vector 
magnitude  combinations  of  data  from  Ringdal  and  Fyen  (1991). 


statistic  greater  than  or  equal  to  actual  values  computed,  Tq,  for  events  24,  25,  26,  and 
27,  and  for  all  15  possible  magnitude  vector  combinations.  These  four  events  were 
the  only  ones  with  probabilities  less  than  0.05.  Recall  that  T*  is  a  measure  of  the 
residual  of  an  event  relative  to  the  sample  random  error  variance;  a  small  probability 
represents  a  large  residual,  discordant  with  the  rest  of  the  data. 

Event  24  was  detected  as  an  outlier  for  case  5,  i.e.,  for  the  case  with  m  = 
(mi, m2).  For  this  ceise,  event  24  is  the  only  small  event  since  m2  v/as  not  observed. 
When  event  24  is  also  removed  from  the  calibration  data,  to  be  tested  as  an  outlier, 
the  lines  of  best  fit  are  given  by  the  dotted  lines  in  Figure  2.  Thus,  with  no  other 
small  events  to  influence  the  slope  and  intercept  estimates,  it  is  not  surprising  that 
event  24  was  detected  as  an  outlier.  In  fact,  event  24  was  very  nearly  identified  as  an 
outlier  for  cases  11  and  12,  which  also  involve  mi  and  m2. 

Event  25  was  detected  as  on  outlier  for  case  5  because  the  residuals  contra¬ 
dict  the  calibration  results.  The  mean  magnitudes  computed  from  the  calibration  data 
are  mi  =  5.710  and  m2  =  5.690,  the  random  error  standard  deviation  estimates  are 


20 


0.072  and  0.079,  and  the  random  error  correlation  estimate  is  very  high  (0.876).  The 
residuals  for  event  25,  however,  are  -0.064  and  +0.065,  respectively,  which  is  almost  a 
two  standard  deviation  difference.  Furthermore,  the  95%  log  yield  confidence  interval 
for  this  case  is  undefined,  i.e.,  there  are  no  real  values  of  w  for  which  the  inequality 
in  (18)  is  satisfied. 

Event  27  was  identified  as  an  outlier  for  cases  6,  7,  and  9.  Cases  7  and 
9  involve  (world-wide  mj)  for  which  the  magnitude  was  6.04,  the  largest  of  all 
magnitudes  mezisured  for  the  18  events.  Likewise,  the  NAO  P  coda  magnitude  for  this 
event  was  5.952,  the  second  largest  observation  for  all  magnitude  types.  Similar  to 
event  25  for  case  5,  the  mean  calibration  magnitudes  are  very  close  and  the  random 
error  correlation  estimates  are  high,  but  the  residuals  are  quite  different.  Thus,  the 
outliers  detected  for  event  27  are  the  result  of  large  residuals  for  m3  and  m^,  and  the 
contradiction  between  the  calibration  data  and  this  event  for  the  different  magnitude 
types.  Note  that  for  event  27,  c<ises  6  and  7,  there  is  no  real  value  of  w  for  which  the 
inequality  in  (18)  is  satisfied. 

Event  26  was  determined  to  be  an  outlier  for  all  cases  not  involving 
(GRF  Lg).  Note  that  mj  was  not  observed  for  event  26.  Apart  from  cases  5,  6,  7, 
and  9,  event  26  was  the  only  outlier  detected.  In  fact,  the  probabilities  that  this 
event  is  consistent  with  the  rest  of  the  data  are  noticeably  smaller  than  for  any  other 
event.  This  strongly  suggests  that  event  26  was  unusual.  There  are  a  number  of 
possible  explanations.  Event  26  might  have  been  decoupled  explosion.  In  fact,  a 
relatively  small  decoupling  factor  would  account  for  the  large  residuals  and  the  high 
slope  estimates  when  event  26  was  included  in  the  calibration  (solid  lines.  Figure  2). 
Complicated  near-source  topography  and  propagation  effects  could  also  be  responsible. 
The  actual  cause  of  the  large  residuals  for  event  26  warrants  further  study. 

Since  event  26  has  been  shown  to  be  an  outlier  (for  whatever  reason),  we 
suggest  that  the  best  calibration  estimates,  based  on  the  magnitudes  from  Ringdal 
and  Fyen  (1991),  are  those  computed  with  event  26  excluded,  i.e.,  the  intercept,  slope 
and  standard  deviation  estimates  corresponding  to  the  dashed  lines  in  Figure  2  and 
the  correlation  coefficient  estimates  in  Table  5. 


4.3  Confidence  Intervals  and  F-Numbers. 


The  100(1  —  a)%  confidence  interval,  presented  in  Section  3.2,  may  be  used 
to  estimate  future  yields,  and  to  test  compliance  to  the  TTBT.  Unfortunately,  a 
single  unified  yield  estimator  whose  distribution  is  known  could  not  be  obtained  while 
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treating  the  slope  parameters  as  unknown.  Clearly,  the  large  uncertainty  in  the  slope 
estimates  warrants  this  treatment.  Hence,  we  have  a  confidence  interval  that  may  be 
applied  to  15  distinct  multivariate  magnitude  combinations.  Also,  since  the  calibration 
parameters  are  unknown,  we  presented  a  definition  of  the  F-number  which  is  a  random 
variable.  There  are  two  questions  we  hope  to  answer  here: 

•  Which  multivariate  magnitude  combination  provides  the  “best”  confidence  in¬ 
terval,  based  on  the  available  calibration  data? 


•  What  are  the  means  and  standard  deviations  of  the  F-numbers? 

To  address  these  questions  we  applied  jackknifing  to  the  data.  The  jackknife 
is  a  resampling  scheme  in  which  each  event  is  successively  removed  from  the  calibration 
data  and  used  as  a  “new”  event.  Hence,  a  confidence  interval  and  random  F-number 
sample  may  be  obtained  for  each  event.  Sample  means  and  standard  deviations  may 
then  be  computed.  In  the  following  analysis,  we  set  a  —  0.05. 

If  jackknifing  is  routinely  applied  to  all  of  the  data,  several  problems  occur. 
First,  the  F-numbers,  treating  event  24  as  the  new  one,  are  ill-defined  for  all  cases 
involving  m2,  since  neither  of  the  small  events,  24  and  26,  were  used  in  the  calibration. 
For  cases  11,  12,  14,  and  15,  the  confidence  intervals  are  open  and,  hence,  the  sample 
F-numbers  are  infinite.  For  ceises  2,  5,  8,  and  9,  the  confidence  intervals  are  closed, 
but  the  sample  F-numbers  are  unusually  large  (>  200).  Second,  the  95%  confidence 
interval  and  F-number  for  event  25  does  not  exist  for  case  5,  nor  do  they  exist  for 
event  27  for  cases  6  and  7;  the  roots  of  the  quadratic  inequality  are  complex.  These 
events  were  identified  as  outliers  for  these  cases.  Third,  since  event  26  is  an  outlier 
for  every  case,  none  of  the  95%  confidence  intervals  include  the  actual  yield.  In  fact, 
a  99%  confidence  interval  would  include  the  actual  yield  of  event  26  in  only  one  C2ise 
(case  13).  For  these  reasons,  the  results  presented  here  were  obtained  by  omitting 
event  26  from  the  analysis  and  without  jackknifing  on  event  24,  i.e.,  event  24  was 
not  removed  from  the  calibration  data  to  be  used  as  a  new  event.  Also,  cases  with 
undefined  confidence  intervals  were  omitted. 

Table  8  lists  the  F-number  mean  and  standard  deviation  estimates,  denoted 
by  F  and  dp,  for  all  vector  magnitude  combinations.  The  third  column  lists  the  num¬ 
ber  of  “new”  events  used  to  obtain  the  jackknife  estimates.  With  event  26  excluded 
from  the  calibration  data,  the  95%  confidence  intervals  are  undefined  (complex)  for 
event  25,  case  5,  and  event  27,  case  7.  The  confidence  intervals  exist,  but  do  not 
contain  the  actual  yield  for  event  27,  cases  4,  6,  and  9  (and  very  nearly  for  case  3). 
These  cases  involve  m3  and  17x4  for  which  the  magnitudes  measured  for  this  event  were 
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Table  8.  Jackknife  estimates  of  the  F-nurhber  means  and  standard  deviations 
for  all  vector  magnitude  combinations  of  data  from  Ringdal  and  Fyen 
(1991). 


Case 

Magnitude  Vector 

^  Samples 

F 

Op 

1 

mi 

15 

1.543 

0.026 

2 

m2 

16 

1.578 

0.034 

3 

ms 

15 

1.637 

0.038 

4 

m4 

16 

1.815 

0.041 

5 

(mi,  m2) 

14 

1.708 

0.058 

6 

(mi,  m3) 

15 

1.676 

0.176 

7 

(mi,m4) 

14 

1.735 

0.072 

8 

(m2,  m3) 

15 

1.626 

0.098 

9 

(m2,m4) 

16 

1.684 

0.122 

10 

(m3,m4) 

15 

1.820 

0.111 

11 

(mi,  m2,  m3) 

15 

1.788 

0.170 

12 

(mi,  m2,  mi) 

15 

1.799 

0.195 

13 

(mi,  ms,  mi) 

15 

1.880 

0.146 

14 

(m2,  ms,  mi) 

15 

1.807 

0.092 

15 

(mi,  m2,  m3,  mi) 

15 

2.013 

0.187 

higher  than  those  for  events  with  similar  yields.  Last,  the  95%  confidence  interval  for 
event  34,  case  3,  does  not  include  the  actual  yield.  The  NAO  P  coda  magnitude  for 
this  event  is  the  smallest  of  all  the  events  in  this  range. 

The  results  in  Table  8  demonstrate  the  following; 

•  Confidence  intervals  based  on  NAO  Lg  magnitudes  are  the  shortest,  followed 
closely  by  those  bzised  on  GRF  Lg  magnitudes.  These  magnitudes  were  the  ones 
with  the  least  random  scatter. 

•  There  is  no  advantage  here  in  using  intervals  based  on  multivariate  magnitude 
vectors.  This  may  be  a  result  of  the  high  random  error  correlations.  The  lowest 
mean  F-numbers,  for  a  magnitude  vector  with  p  >  1,  were  obtained  for  cases  6, 
8  and  9  for  which  the  random  error  correlations  are  the  smallest  (Table  5). 

•  The  sample  F-number  standard  deviations  for  the  first  five  cases  are  reasonably 
small,  suggesting  that  the  corresponding  sample  means  are  stable  estimates. 
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SECTION  5 


RESULTS:  ISRAELSSON  MAGNITUDE  DATA 


5.1  Calibration  Results. 


Using  the  Soviet  network  data  in  Table  3,  we  have  computed  calibration 
parameter  estimates.  Figure  4  shows  scatter  plots  of  the  data  and  lines  of  best  fit. 
As  before,  three  lines  of  best  fit  were  computed  using  all  17  events  (solid),  all  except 
event  26  (dashed),  and  all  except  events  24  and  26  (dotted).  Events  24  and  26  are 
represented  by  the  solid  square  and  triangle  markers,  respectively,  and  the  intercept, 
slope  and  standard  deviation  estimates  are  given  in  the  legends  above  each  frame. 
R2mdom  error  correlation  coefficient  estimates  are  given  in  Tables  9-11  for  these  three 
Ccises. 


There  are  several  noteworthy  points  regarding  these  results.  First,  the  mag¬ 
nitudes  based  on  RMS  initial  P  amplitudes  exhibit,  surprisingly,  the  least  random 
scatter  with  a  sample  random  error  standard  deviation  of  0.095,  as  compared  to  0.111, 
0.123,  and  0.121  for  the  magnitudes  bcised  on  the  RMS  amplitudes  of  P  coda,  S  with 
coda,  and  Lg,  respectively  (Figure  4).  Israelsson  (1992)  noted  that  there  were  clearly 
developed  Lg  phases  at  only  two  of  the  stations,  APA  and  ARU,  which  were  two  of 
the  three  closest  to  NNZ.  Poorly  developed  Lg  ph2ises  at  the  other  stations  could  be 
due  to  Lg  blockage  in  the  Barents  and  Kara  Seas  (Baumgardt,  1990;  Israelsson,  1992). 
Israelsson  (1992)  stated,  however,  that  there  was  sufficient  SNR  to  compute  RMS  Lg 
amplitudes  even  for  the  two  smallest  explosions.  The  results  obtained  here  suggest 
that  Lg  blockage  and  low  SNR  (relative  to  that  of  the  initial  P  and  P  coda  phases) 
may  have  some  impact  or,  possibly,  that  near-source  effects  and  surface  interactions 
are  responsible.  Further  study  is  needed  to  understand  this  result. 

Second,  the  slope  estimates,  based  on  all  of  the  events  (solid  lines),  are  again 
larger  than  expected  and,  hence,  the  intercept  estimates  are  smaller  than  expected. 
Note  that  when  event  26  is  omitted,  the  slope  estimates  decrease  significantly  (Figure 
4).  The  sample  random  error  standard  deviations  also  decrease  for  all  four  magnitudes 
when  this  event  is  omitted.  When  both  events  26  and  24  are  omitted,  the  slope 
estimates  for  m4(Scoda)  and  mi(Lg),  in  particular,  are  indeterminate  due  to  the  tight 
clustering  of  the  remaining  events.  Third,  the  correlation  coefficient  estimates  are 
remarkably  high  for  all  cases  (Tables  9-11). 


24 


mb(Scoda) 


N  =  17  a  =  3.383  b  =  1.141  a  =  0.111 

N  =  18  a  =  3.841  b  =  0.821  a  =  0.092 

N  =  15  a  =  4.778  b  =  0.471  a  =  0.091 


N  =  17  a  =  3.248  b  =  1.200  a  =  0.121 

N  =  16  a  =  3.787  b  =  0.942  a  =  0.095 

N  =  15  a  =  5.927  b  =  -.088  a  =  0.075 


w  (log  kt)  w  (log  kt) 


Figure  4.  Plots  of  Soviet  network  magnitudes,  based  on  RMS  amplitutes  of 
initial  P,  P  coda,  S  with  coda,  and  Lg,  versus  the  yields  of  17  nuclear 
tests  at  NNZ.  The  solid,  dashed  and  dotted  lines  in  each  frame 
represent  the  lines  of  best  fit  using  all  of  the  available  17  events, 
all  except  event  26  (solid  triangle),  and  all  except  events  26  and  24 
(solid  square).  The  parameter  estimates  are  given  in  the  legends 
above  the  plots. 


N  == 

17 

a 

=  3.241 

b  = 

1.201 

o  = 

0.123 

N 

18 

a 

=  3  851 

b  = 

0.908 

a  = 

0.087 

IS 

a 

=  5.387 

b  = 

0.179 

o  = 

0.078 

N  ■ 

17 

a 

=  3.278 

b  » 

1.197 

o  ~ 

0.095 

N  = 

16 

a 

=  3.745 

b  = 

0.972 

0  ~ 

0.068 

N 

15 

a 

-  4.018 

b  = 

0.842 

o  ~ 

0.070 

25 


Table  9.  Random  error  correlation  coefficient  estimates  for  Soviet  network 
magnitude-log  yield  relations,  based  on  the  magnitude  data  from 
Israelsson  (1992). 


Magnitude 

mifP) 

mfcfPcoda) 

mjfScoda) 

mtfLg) 

m4(P) 

1.000 

0.921 

0.877 

0.870 

m4(Pcoda) 

1.000 

0.916 

0.891 

m4(Scoda) 

1.000 

0.953 

y»>(Lg) _ 

1.000 

Table  10.  Correlation  coefficient  estimates  excluding  event  26. 

Magnitude 

mfcfP) 

mtCPcoda) 

mtfScoda) 

mfcfLg) 

mj(P) 

1.000 

0.881 

0.740 

0.761 

m4(Pcoda) 

1.000 

0.873 

0.823 

m4(Scoda) 

1.000 

0.920 

^t(Lg) _ 

1.000 

Table  11.  Correlation  coefficient  estimates  excluding  event  24  and  26. 


Magnitude 

m*(P) 

mjfPcoda) 

mjfScoda) 

»n6(P) 

1.000 

0.892 

0.794 

0.906 

m4(Pcoda) 

1.000 

0.878 

0.869 

m4(Scoda) 

1.000 

0.903 

1.000 
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Table  12.  Probability  of  T*  >  T*  for  events  detected  as  outliers  for  all  vector 
magnitude  combinations  of  data  from  Israelsson  (1992). 


Case 

Magnitude  Vector 

Event 

24 

26 

38 

ms 

0.010 

me 

my 

0.038 

0.001 

mg 

0.006 

5 

(ms,  me) 

0.012 

0.006 

6 

(ms,  my) 

0.038 

0.005 

(ms,  mg) 

0.013 

0.008 

(me,  my) 

0.005 

(me,  mg) 

0.027 

10 

(my,  mg) 

0.037 

0.005 

0.016 

11 

(ms,  me,  my) 

0.003 

12 

(ms,  me,  mg) 

0.015 

13 

(ms,  my,  mg) 

0.011 

0.010 

0.047 

14 

(me,  my,  mg) 

0.014 

0.047 

15 

(ms,  me,  my,  mg) 

0.005 

0.005 

5.2  Outlier  Analysis  Results. 

We  have  applied  the  outlier  test,  described  in  Section  3.3,  to  these  data. 
Table  12  lists  the  probabilities  of  obtaining  values  of  greater  than  or  equal  to  the 
actual  values  computed,  Tq,  for  events  24,  26,  and  38,  and  for  all  distinct  magnitude 
vector  combinations.  These  three  events  were  the  only  ones  with  probabilities  less 
than  0.05.  The  magnitudes  based  on  the  RMS  amplitudes  of  initial  P,  P  coda,  S  with 
coda,  and  Lg  are  denoted  by  ms,  m^,  mj,  and  mg,  respectively. 

As  in  Section  4.2,  event  26  is  an  outlier  for  all  cases  because  of  the  large 
residuals.  Event  24  is  also  an  outlier  for  10  of  the  15  cases.  This  is  caused  by  the  fact 
that  the  calibration  estimates  are  very  poor  when  event  26  is  the  only  small  event. 
This  also  occurred  for  the  analysis  based  on  the  magnitudes  from  Ringdal  and  Fyen 
(1991),  but  not  to  this  degree  since  magnitudes  for  event  42  were  measured.  Event 
42  has  the  third  smallest  yield  and  helped  to  establish  reasonable  slope  and  intercept 
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estimates.  Soviet  network  magnitudes  were  not  available  for  this  event.  If  other  small 
events  were  present  in  the  data  set,  it  is  doubtful  that  event  24  would  be  an  outlier. 

Event  38  is  an  outlier  for  cases  10,  13,  and  14,  which  all  involve  magnitudes 
m^  and  mg.  The  mean  magnitudes  computed  from  the  calibration  data  are  my  = 
mg  =  5.632,  the  random  error  sample  standard  deviations  are  0.127  and  0.121,  and 
the  random  error  correlation  estimate  is  0.973.  The  residuals  for  event  38  are  0.032 
and  0.131,  respectively.  Note  that  for  event  38,  case  10,  there  is  no  real  value  of  w  for 
which  the  inequality  in  (18)  is  satisfied. 

The  preceeding  analysis  indicates  that  event  26  is  also  unusual  for  this  data 
set.  Hence,  we  suggest  that  the  best  calibration  estimates,  based  on  the  magnitudes 
from  Israelsson  (1992),  are  those  computed  with  event  26  excluded,  i.e.,  the  intercept, 
slope  and  standard  deviation  estimates  corresponding  to  the  dashed  lines  in  Figure  4 
and  the  correlation  coefficient  estimates  in  Table  10. 


5.3  Confidence  Intervals  and  F-Numbers. 


Here  we  present  the  jackknife  estimates  of  the  F-number  means  and  standard 
deviations,  based  on  the  magnitude  data  from  Israelsson  (1992).  As  before,  a  =  0.05, 
event  26  was  excluded  from  the  analysis,  and  the  jackknife  was  not  applied  to  eve^ 
24.  Table  13  lists  the  sample  F-number  means  and  standard  deviations,  denoted  by  F 
and  djp,  for  all  vector  magnitude  combinations.  The  third  column  lists  the  number  of 
samples  used  in  the  jackknife  procedure.  The  95%  confidence  intervals  au^e  undefined 
for  event  28,  cases  10  and  13.  The  confidence  intervals  exist,  but  do  not  contain  the 
actual  yields  for  event  38,  case  14,  event  28,  cases  1,  2,  and  8,  and  event  34,  case  1. 

The  results  in  Table  13  demonstrate  the  following: 

•  Confidence  intervals  based  on  initial  P  magnitudes,  mg,  are  the  shortest,  followed 
by  those  based  on  the  p  =  2  magnitude  vectors  (mg,  mg),  (mg,  my)  and  (mg,  mg). 
The  magnitudes  based  on  initial  P  RMS  amplitudes  exhibit  the  least  random 
scatter.  The  sample  mean  F-number  for  mg  is  also  smaller  than  those  based  on 
the  data  from  Ringdal  and  Fyen  (1991).  (Compare  Tables  8  and  13.) 

•  There  is  little  advantage  in  using  intervals  based  on  multivariate  magnitude  vec¬ 
tors.  However,  since  the  sample  random  error  standard  deviation  is  significantly 
smaller  for  mg  than  for  the  other  magnitudes,  the  mean  F-numbers  for  the  three 
p  =  2  magnitude  vectors  listed  in  item  1  are  smaller  than  those  for  mg,  my,  and 
mg  individually  even  though  the  random  error  correlation  estimates  are  high. 
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(Cf.  the  parameter  estimates  corresponding  to  the  dcished  lines  in  Figure  4  and 
in  Table  10.) 

•  Most  of  the  sample  F-number  standard  deviations  are  relatively  small,  suggesting 
that  the  samples  means  are  stable  estimates. 


Table  13.  Jackknife  estimates  of  the  F-number  means  and  standard  deviations 
for  all  vector  magnitude  combinations  of  data  from  Israelsson  Fyen 
(1992). 


Case 

Magnitude  Vector 

^  Samples 

F 

Op 

1 

ms 

15 

1.448 

2 

me 

15 

1.720 

3 

mj 

15 

1.672 

4 

ms 

15 

1.726 

0.042 

5 

(ms,  me) 

15 

1.547 

0.089 

6 

(m5,m7) 

15 

1.601 

0.078 

7 

(ms,  ms) 

15 

1.602 

0.062 

8 

(m6,m7) 

15 

1.926 

9 

(me,  ms) 

15 

1.956 

10 

(m7,m8) 

15 

1.984 

BSiH 

11 

(ms,m6,  m7) 

15 

1.647 

12 

(ms,  me,  mg) 

15 

1.697 

13 

(ms,m7,  mg) 

15 

1.780 

H^B 

14 

(m6,m7,  mg) 

15 

2.205 

15 

(ms,  me,  m7,  mg) 

15 

1.773 
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SECTION  6 


CONCLUSIONS 


The  Soviet  yield  data  set  presented  here  provides  the  first  published  set  of 
actual  yields  for  undergroxind  nuclear  tests  conducted  at  Novaya  Zemlya.  We  have 
also  provided  the  first  direct  estimates  of  the  calibration  parameters  for  several  seismic 
magnitudes,  NORSAR  Lg  and  P  coda,  Grafenberg  Lg,  a  world-wide  m^,  and  four 
Soviet  network  magnitudes,  based  on  RMS  P,  P  coda,  S  with  coda,  and  Lg  amplitudes. 
As  noted  in  Sections  1  and  2.1,  the  accuracy  of  the  yield  data  must  be  resolved  before 
a  high  level  of  confidence  should  be  placed  in  these  results. 

The  calibration  results  in  Sections  4.1  and  5.1  showed  that  the  slope  and 
intercept  estimates  depend  highly  on  the  two  smallest  events.  In  Sections  4.2  and  5.2, 
we  showed  that  the  residuals  for  one  of  these  tests,  event  26,  were  outliers  for  every 
relevant  magnitude  combination.  In  addition,  event  26  was  responsible  for  yielding 
unusually  high  slope  and  low  intercept  estimates.  For  these  reasons,  our  best  estimates 
of  the  calibration  parameters  are  those  in  Tables  14-17,  which  were  computed  with 
event  26  excluded. 

Treating  all  of  the  calibration  parameters  as  unknown,  we  provided  a  classi¬ 
cal  95%  confidence  interval,  by  which  future  yields  may  be  estimated  and  compliance 
to  a  treaty  threshold  may  be  tested.  An  F-number  was  defined  in  terms  of  the  interval, 
such  that  (y /F,  y  X  F)  hzis  a  95%  probability  of  including  the  actual  yield,  where  y 
is  the  central  value  of  the  interval.  This  definition  was  useful  here  because  we  do  not 
have  a  log  yield  estimator  whose  distribution  is  known.  Sample  means  and  standard 
deviations  of  the  F-numbers  were  computed  for  various  magnitude  combinations  by 
jackknifing.  The  means  ranged  from  1.448  to  2.205  (Tables  8  and  13),  depending  on 
the  magnitude  combination. 

Since  event  26  wzis  shown  to  be  an  outlier,  there  was  only  one  small  event 
(<  30  KT)  used  to  estimate  the  slopes  and  intercepts;  the  remainder  were  clustered 
between  97  and  139  KT.  Thus,  the  accuracy  of  the  slope  and  intercept  estimates 
depend  largely  on  event  24.  A  Bayesian  approach,  that  can  incorporate  a  priori 
expert  opinion  with  the  magnitude-yield  data,  may  prove  to  be  the  best  method  for 
estimating  the  yield  of  a  new  event  and  testing  compliance.  It  provides  a  statistical 
basis  for  incorporating  a  priori  information  for  cases  in  which  the  data  are  insufficient 
to  compute  reliable  calibration  estimates. 
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Shumway  and  Der  (1990)  have  developed  such  an  approach  by  which  the 
confidence  intervals  may  be  computed  analytically.  Their  analytic  approach  is  some¬ 
what  limited,  however,  by  the  form  of  the  Bayesian  prior  distribution  for  the  unknown 
random  error  covariance  matrix.  Fisk  et  al.  (1992)  have  also  developed  a  Bayesian 
method  which,  in  addition  to  the  other  information,  can  incorporate  information  avail¬ 
able  from  no-yield  magnitude  data  (irrelevant  for  NNZ  where  we  actually  have  yield 
data  for  more  events  than  we  have  magnitude  data).  In  addition,  the  approach  allows 
arbitrary  specification  of  the  Bayesian  prior  distribution  for  the  covariance  matrix. 
Yield  estimates  and  the  critical  value  of  the  hypothesis  test  of  TTBT  compliance 
are  computed  numerically.  Gray  et  al.  (1992)  have  extended  the  method  of  Fisk  et 
al.  (1992)  to  treat  unknown  slope  parameters. 

With  the  possibility  of  resumed  testing  at  NNZ  in  late  1992,  a  simulation 
to  assess  the  performances  of  the  various  yield  estimation  methods  warrants  future 
investigation. 
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Table  14.  Best  intercept,  slope,  and  random  error  standard  deviations 
estimates  for  NAO  Lg,  GRF  Lg,  NAO  P  coda  and  world-wide  mj 
magnitudes. 


Magnitude 

a 

b 

a 

NAO  Lg 

3.975 

0.865 

0.071 

GRF  Lg 

3.800 

0.942 

0.081 

NAO  Pcoda 

3.672 

1.012 

0.093 

mj 

3.431 

1.126 

0.124 

Table  15.  Best  random  error  correlation  coefficient  estimates  for  NAO  Lg,  GRF 
Lg,  NAO  P  coda  and  world-wide  mj  magnitudes. 


Magnitude 

NAO  Lg 

GRF  Lg 

NAO  Pcoda 

mi, 

NAO  Lg 

1.000 

0.787 

0.664 

0.557 

GRF  Lg 

1.000 

0.441 

0.383 

NAO  Pcoda 

1.000 

0.946 

mj 

1.000 
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Table  16.  Best  intercept,  slope  and  random  error  standard  deviation  estimates 
for  Soviet  network  magnitudes  based  on  initial  P,  P  coda,  S  with 
coda,  and  Lg  RMS  amplitudes. 


Magnitude 

a 

b 

d 

m6(P) 

3.745 

0.972 

0.068 

mi,(  Pcoda) 

3.841 

0.921 

0.092 

mj,(Scoda) 

3.851 

0.908 

0.087 

m6(Lg) 

3.787 

0.942 

0.095 

Table  17.  Best  random  error  correlation  coefficient  estimates  for  Soviet 
network  magnitudes  based  on  initial  P,  P  coda,  S  with  coda,  and 
Lg  RMS  amplitudes. 


Magnitude 

m6(P) 

m6  (Pcoda) 

mj(Scoda) 

^«»(Lg) 

1.000 

0.881 

0.740 

0.761 

mt(Pcoda) 

1.000 

0.873 

0.823 

m6(Scoda) 

1.000 

0.920 

mt{Lg) 

1.000 
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